package NumTheory;

import java.math.BigDecimal;
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.Map;
import java.util.Random;

public class NumTheory implements NumTheoryInterface {
    @Override
    // 算法输入为正整数n，输出为n以内的所有素数
    public BigInteger[] eratosthenes(BigInteger n) {
        // 由于是BigInteger，会超标，所以用HashMap代替数组
        Map<BigInteger, Boolean> map = new HashMap<>();
        // 每个元素最初的默认值为True
        for (BigInteger i = BigInteger.ONE; i.compareTo(n) <= 0; i = i.add(BigInteger.ONE)) {
            map.put(i, Boolean.TRUE);
        }
        // 计算根号n的值m
        BigInteger m = BigNumSqrt(n);
        for (BigInteger i = BigInteger.valueOf(2); i.compareTo(m) <= 0;
             i = i.add(BigInteger.ONE)) {
            if (map.get(i)) {
                // 将map中i的倍数下标标记为false（但是不包括i的1倍）
                for (BigInteger j = new BigInteger(String.valueOf(i)).multiply(BigInteger.valueOf(2)); j.compareTo(n) <= 0;
                     j = j.add(i)) {
                    map.put(j, Boolean.FALSE);
                }
            }
        }
        // 返回所有标记为true的元素
        ArrayList<BigInteger> arrayList = new ArrayList<>();
        for (BigInteger i = BigInteger.valueOf(2); i.compareTo(n) <= 0;
             i = i.add(BigInteger.ONE)) {
            if (map.get(i)) {
                arrayList.add(i);
            }
        }
        return arrayList.toArray(new BigInteger[0]);
    }

    @Override
    public BigInteger euclid(BigInteger a, BigInteger b) {
        // 只能处理两个正整数a和b
        if (a == null || b == null || a.equals(BigInteger.ZERO) || b.equals(BigInteger.ZERO)) {
            return BigInteger.ZERO;
        }
        while (b.compareTo(BigInteger.ZERO) != 0) {
            BigInteger mod = a.mod(b);
            a = b;
            b = mod;
        }

        return a;
    }

    @Override
    public BigInteger[] exEuclid(BigInteger a, BigInteger b) {
        // 只能处理两个正整数a和b
        if (b.equals(BigInteger.ZERO)) {
            return new BigInteger[]{a, BigInteger.ONE, BigInteger.ZERO};
        } else {
            BigInteger[] res = exEuclid(b, a.mod(b));
            BigInteger x = res[2];
            BigInteger y = res[1].subtract(a.divide(b).multiply(res[2]));
            return new BigInteger[]{res[0], x, y};
        }
    }

    @Override
    public BigInteger fastPow(BigInteger base, BigInteger expo, BigInteger p) {
        BigInteger result = BigInteger.ONE;
        while (expo.compareTo(BigInteger.ZERO) > 0) {
            if (expo.mod(BigInteger.valueOf(2)).intValue() == 1) {
                result = result.multiply(base).mod(p);
            }
            expo = expo.divide(BigInteger.valueOf(2));
            base = base.multiply(base).mod(p);
        }

        return result;
    }

    @Override
    public BigInteger crt(BigInteger[] a, BigInteger[] p, int num) {
        // 余数数组是a，模数数组是p，方程个数是num
        BigInteger P = BigInteger.ONE;
        for (int i = 0; i < num; i++) {
            P = P.multiply(p[i]);
        }
        BigInteger result = BigInteger.ZERO;
        // 生成逆元数组
        BigInteger[] PInv = new BigInteger[num];
        for (int i = 0; i < num; i++) {
            PInv[i] = P.divide(p[i]).modInverse(p[i]);
            result = result.add(P.divide(p[i]).multiply(PInv[i]).multiply(a[i]));
        }

        return result.mod(P);
    }

    @Override
    public boolean MRtest(BigInteger n) {
        // 找到满足 n-1 = 2^k × q
        BigInteger k = BigInteger.ZERO;
        BigInteger q;
        BigInteger nsub1 = n.subtract(BigInteger.ONE);
        while (nsub1.mod(new BigInteger("2")).compareTo(BigInteger.ZERO) == 0) {
            nsub1 = nsub1.divide(new BigInteger("2"));
            k = k.add(BigInteger.ONE);
        }
        q = nsub1;

        // 在(1,n)中找到一个随机数 a
        Random random = new Random();
        BigInteger a;
        do {
            a = new BigInteger(n.bitLength(), random);
        } while (a.compareTo(n) >= 0 && a.compareTo(BigInteger.ONE) <= 0);

        if (fastPow(a, q, n).compareTo(BigInteger.ONE) == 0) return true;

        BigInteger twoPowI = BigInteger.ONE;
        for (BigInteger i = BigInteger.ZERO; i.compareTo(k) < 0; i = i.add(BigInteger.ONE)) {
            if (fastPow(a, twoPowI.multiply(q), n).compareTo(n.subtract(BigInteger.ONE)) == 0) return true;
            twoPowI = twoPowI.multiply(BigInteger.valueOf(2));
        }

        return false;
    }

    /**
     * Miller-Rabin素性检测算法，但是增加了轮数，用以提高准确度
     *
     * @param n    待检测的大整数
     * @param time 检测的轮数
     * @return 返回是否是素数
     */
    public boolean MRtest(BigInteger n, int time) {
        boolean flag = true;
        while (time-- >= 1) {
            flag = MRtest(n);
            if (!flag) {
                break;
            }
        }
        return flag;
    }


    // 大数开方，使用牛顿迭代法近似逼近
    // 如果看不懂，可以看看博客
    // https://blog.csdn.net/qq_38701868/article/details/109145801
    public BigInteger BigNumSqrt(BigInteger BigNum) {
        if (BigNum == null || BigNum.compareTo(BigInteger.ZERO) == 0) {
            return BigInteger.ZERO;
        }
        // 求解数N的平方根x，也就是x^2=N，也就是求解公式 y=x^2-N 的零解
        // 1.取一点x_i，那么曲线上的点为(x_i, x_i^2-N)，切线斜率是2x_i
        // 2.此时切线为 y-(x_i^2-N)=2x_i(x-x_i)
        // 3.整理切线得 y=x_i^2-N+2x_i(x-x_i)=x_i^2-N+2x_ix-2x_i^2=2x_ix-x_i^2-N
        // 4.与x轴交点为((x_i^2+N)/(2x_i), 0) = (0.5(x_i+N/x_i), 0)
        // 5 x_{i+1} = 0.5(x_i+N/x_i)
        BigDecimal newBigNum = new BigDecimal(BigNum);
        BigDecimal xi = BigDecimal.ONE;
        // 如果差距大于0.001，那还要继续逼近
        while ((xi.multiply(xi).subtract(newBigNum)).abs().compareTo(BigDecimal.valueOf(0.001)) > 0) {
            // i++;
            BigDecimal s = xi.add(newBigNum.divide(xi, 2000, BigDecimal.ROUND_HALF_UP));
            xi = s.divide(BigDecimal.valueOf(2), 2000, BigDecimal.ROUND_HALF_UP);
        }
        return new BigInteger(xi.toString().split("\\.")[0]);
    }
}
